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Abstract 


In  this  paper  /pa  describe'pDRPACK,  a  software  package  for  the  Or¬ 
thogonal  Distance  R^gr^jrrq  (OPR)  problem.  This  software  implements 
the  algorithm  ftesnibcd  in  (BogBSflfi)  for  finding  the  set  of  parameters  that 
minimise  the  sum  of  the  squared  orthogonal  distances  from  a  set  of  ob¬ 
servations  to  a  curve  determined  by  the  parameters.  It  can  also  be  used 
to  solve  the  ordinary  nonlinear  least  squares  problem.  The  ODR  proce¬ 
dure  has  application  to  curve  fitting,  and  to  the  errors  in  variables  problem 
in  statistics.  The  algorithm  implemented  is  an  efficient  and  stable  trust- 
region  ( Levenberg-Mar quardt)  procedure  that  exploits  the  structure  of  the 
problem  so  that  the  computational  cost  per  iteration  is  equal  to  that  for 
the  same  type  of  algorithm  applied  to  the  ordinary  nonlinear  least  squares 
problem.  The  package  allows  a  general  weighting  scheme,  provides  for  fi¬ 
nite  difference  derivatives,  and  contains  extensive  error  checking  and  report 

generating  facilities.  -  . . . 
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1  Introduction 


In  [BogBS&5]  the  authors  provide  an  efficient  and  stable  algorithm  for  the 
orthogonal  distance  regression  (ODR)  problem.  This  problem  is  defined  as 
follows.  Let  (x*,g!(),  t  =  be  a  given  set  of  data  where  €  Rm 

and  pi  €  I1.  Suppose  that  the  values  of  m  ana  function  of  x,  and  a  set 
of  unknown  parameters  fi  €  RF,  but  that  both  the  ft  and  the  x<  contain 
actual,  but  unknown,  errors  ef  eRl  and  6*  €  R",  respectively.  Then  the 
observed  value  of  y,  satisfies 

*  = /(*  +  $<•; /**)-«?  s*l,...,» 

for  some  actual  but  again  unknown  value  0*.  (Note  that  we  have  chosen 
the  sign  of  a  to  be  negative  for  convenience.)  The  orthogonal  distance 
regression  problem  is  to  approximate  by  finding  the  fi  for  which  the 
sum  of  the  squares  of  the  n  orthogonal  distances  from  the  curve  /(x;0)  to 
the  n  data  points  is  minimised.  This  ie  accomplished  by  the  minimisation 
problem 

***  im 1 

subject  to  the  constraints 

*  =  /(*  +  />)  -  ti  »  =I,...,n. 

Since  these  constraints  are  linear  in  e,,  we  can  eliminate  them  and  the  t,, 
thereby  obtaining 


<5j»E  [(/(*  +  *;«-*)’  +  *7*1  •  (i) 

We  obtain  the  final  form  of  the  problem  to  be  solved  by  allowing  a 
general  weighting  scheme.  Let  w,,  i*  1, ...  ,n,  be  a  set  of  non-negative 
numbers  and  let  Di  €  1**",  »  =  l, ...  ,n,  be  positive  diagonal  matrices. 
Then  (l)  is  generalised  to  the  weighted  ortkofonml  distance  regression  prob¬ 
lem  n 

■WE"?  ((/(*  +  *;«-»)’  +  (J) 
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This  weighting  scheme  enables  ns  to  apply  our  procedure  to  problems  where 
U  and  St  hare  different  variances,  and  to  situations  where  various  observe- 
tions  should  be  weighted  differently. 

The  algorithm  developed  in  [BogBSSft]  is  a  trust-region,  Levenberg- 
Marquardt  method  that  exploits  the  structure  of  the  problem  to  obtain  a 
procedure  that  is  both  stable  and  efficient.  In  the  case  of  ordinary  least 
squares  (OLS),  i.e.,  where  all  of  the  St  are  assumed  to  be  sero,  the  trust 
region,  Levenberg-Marquardt  method  requires  0(np*)  operations  per  it¬ 
eration.  (See,  e.g.,  [DenSSS]  and  [Mor77].)  The  order  of  operations  per 
iteration,  and  the  constant  for  the  highest  order  term,  is  the  same  for  the 
algorithm  developed  in  (BogBSSft]  as  for  OLS.  Note  that  a  straight  forward 
use  of  an  OLS  algorithm  on  (2)  would  require  0(n(n  +  p)2)  operations  per 
iteration  (see  [BogBSSft])  which  is  dearly  prohibitive  for  large  values  of  n. 

ODRPACK  implements  the  algorithm  presented  in  |BogB8Sft].  In  f  2 
we  describe  ODRPACK  in  detail  and  in  $  3  we  discuss  the  installation  of 
the  package.  Appendices  A  and  B  give  a  sample  program  and  its  output, 
respectively.  The  complete  user’s  reference  guide  and  installation  manual 
are  distributed  with  the  package. 

2  Package  Overview 

2.1  Highlights 

ODRPACK  is  a  portable  collection  of  ANSI  77  Fortran  subprograms  for 
fitting  a  model  to  data.  It  is  designed  primarily  for  instances  when  the  in¬ 
dependent  as  well  as  the  dependent  variables  have  significant  errors.  ODR¬ 
PACK  embodies  a  highly  eflkient  algorithm  for  solving  the  weighted  or¬ 
thogonal  distance  regression  problem.  In  addition,  it  can  be  used  to  solve 
the  ordinary  least  squares  problem  where  all  of  the  errors  are  attributed  to 
the  observations  of  the  dependent  variable,  y,. 

ODRPACK  is  designed  to  accommodate  many  levels  of  user  sophistica¬ 
tion  and  problem  difficulty. 

e  It  is  easy  to  use,  providing  two  levels  of  user-control  of  the  com¬ 
putations,  extensive  error  handling  facilities,  comprehensive  printed 


reporta,  complete  online  documentation  and  no  site  restriction*  other 
than  effective  machine  sise. 

•  The  necessary  derivatives  (Jacobian  matrices)  are  approximated  nu¬ 
merically  if  they  are  not  supplied  by  the  user.  In  addition,  the  cor¬ 
rectness  of  user-supplied  derivatives  can  be  verified  by  the  derivative 
checking  procedure  provided. 

•  Both  weighted  and  unweighted  analysis  can  be  performed. 

•  Subsets  of  the  unknowns  can  be  treated  as  constants  with  their  values 
held  fixed  at  their  input  values,  allowing  the  user  to  examine  models 
with  fewer  parameters  without  rewriting  the  model  subprogram. 

•  ODRPACK  has  a  default  scaling  algorithm  that  attempts  to  auto¬ 
matically  accommodate  poorly  scaled  problems  in  which  the  model 
parameters  and/or  unknown  errors  in  the  independent  variables  vary 
widely  in  magnitude.  Alternatively,  the  user  can  supply  appropriate 
scaling  values. 

•  ODRPACK  is  portable  and  is  easily  used  with  other  Fortran  subpro¬ 
gram  libraries.  (See  §  3.3.) 

The  following  sections  discuss  the  use  of  ODRPACK  and  provide  a  brief 
description  of  ODRPACK *s  main  features. 

2.2  Implementation 

3.3.1  Algorithm 

A  full  discussion  of  the  ODRPACK  algorithm  is  found  in  |BogBS*S).  Briefly, 
ODRPACK  implements  a  trust  region  Levenberg-Marquardt  method,  using 
scaling  to  accommodate  problems  in  which  estimated  values  have  widely 
varying  magnitudes.  The  Jacobian  matrices,  i.e.,  the  matrices  of  first  par¬ 
tial  derivatives  of  /(x;/J)  with  respect  to  0  and  x,  are  computed  at  every 
iteration  either  by  finite  differences  or  by  a  user-supplied  subprogram  (see 
§  2.3.3).  The  iterations  are  stopped  when  any  one  of  three  stopping  criteria 
are  met.  Two  of  these  indicate  the  iterations  have  converged  to  a  solution. 
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These  are  evm-of-tqvares  convergence ,  which  indicates  that  the  change  in 
the  sum  of  the  squared  weighted  observational  errors  is  sufficiently  small, 
and  parameter  convergence ,  which  indicates  the  change  in  the  parameters 
is  sufficiently  small.  The  third  stopping  criteria  is  a  limit  on  the  number  of 
iterations. 

2.2.2  Starting  Values 

The  user  must  supply  starting  values  for  the  unknowns  being  estimated,  i.e., 
for  P  and  6.  Users  familiar  with  the  ordinary  nonlinear  least  squares  prob¬ 
lem  are  generally  aware  of  the  importance  of  obtaining  good  starting  values 
for  the  estimated  function  parameters.  It  is  equally  important  here.  Good 
initial  approximations  can  significantly  decrease  the  number  of  iterations 
required  to  find  a  solution;  a  poor  initial  approximation  may  even  prevent 
a  solution  from  being  found  at  all.  Reasonable  initial  approximations  are 
often  available  from  previous  analysis  or  experiments.  When  good  starting 
values  are  not  readily  available,  the  user  may  have  to  do  some  preliminary 
analysis  to  obtain  them  (Him70]. 

Users  who  do  not  provide  scale  information  are  strongly  encouraged  not 
to  use  sero  as  an  initial  approximation  for  any  of  the  function  parameters, 
Pk,  &  =  1, . . .  ,p,  since  a  sero  value  can  result  in  incorrect  scaling  (see  ' 
§  2.2.4).  Setting  the  initial  approximation  to  the  largest  magnitude  which, 
for  the  user’s  problem,  is  effectively  sero,  rather  than  the  actual  value 
sero,  will  help  to  eliminate  scaling  problems  and  possibly  produce  faster 
convergence.  For  example,  if  P\  represents  change  in  a  cost  measured  in 
millions  of  dollars,  then  the  value  10  might  be  considered  “effectively  sero” 
and  an  initial  value  of  /Jf  =  10  is  preferable  to  ft  =  0. 

When  using  orthogonal  distance  regression  it  is  also  important  to  have 
good  starting  values  for  the  estimated  errors,  4,  »  =  1, . . .  ,n.  The  ODR- 
PACK  default  is  to  initialise  6,  to  sero,  which  is  the  most  obvious  initial 
value.  (Note  that  sero  starting  values  for  A  do  not  cause  scaling  problems 
similar  to  those  discussed  above  for  p.)  Initialising  A  to  sero,  however,  is 
equivalent  to  initially  assigning  all  of  the  errors  to  the  dependent  variable 
as  is  done  for  ordinary  least  squares.  While  this  is  quite  adequate  in  many 
cases,  in  others  it  is  not.  A  plot  of  the  observed  data  and  of  the  curve 


described  by  the  model  function  for  the  initial  parameters  may  indicate 
whether  or  not  zero  starting  values  for  Si  are  reasonable  and  may  make  it 
possible  to  visually  determine  better  starting  values  for  some  of  the  For 
example,  if  such  a  plot  shows  that  the  vertical  distance  from  a  data  point 
(z,,y()  to  the  curve  is  far  larger  than  the  orthogonal  distance,  then  Si  should 
probably  not  be  initialized  to  zero.  This  may  occur  near  an  asymptote  or 
near  a  local  minimum  or  maximum.  In  such  cases,  it  is  often  appropriate 
to  initialize  6,  to  the  horizontal  distance  from  the  data  point  to  the  curve. 

2.2.3  Derivative  Handling 

As  was  noted  in  §  2.2.1  the  Jacobian  matrices,  i.e.,  the  matrices  of  first 
partial  derivatives  of  f(x;0)  with  respect  to  each  /?  and  each  component 
of  Xi,  are  required  at  every  iteration.  The  user  may  provide  the  necessary 
derivatives,  or  they  can  be  computed  by  ODRPACK  using  finite  differences. 
Finite  difference  derivatives  generally  cause  very  little  change  in  the  results 
from  those  obtained  using  analytic  derivatives,  provided  sufficient  precision 
is  used  (typically  double  precision  on  a  32-bit  machine). 

Because  coding  errors  are  a  common  problem  with  user-supplied  deriva¬ 
tives,  ODRPACK  has  an  option  to  check  the  validity  of  the  user-supplied 
derivative  code  by  comparing  its  results  to  finite  difference  values  for  the  * 
derivative.  The  derivative  checking  procedure  examines  only  one  row  of  the 
Jacobian  matrices,  and  is  therefore  quite  efficient.  Checking  only  one  row 
is  reasonable  for  regression  models  since  the  same  code  is  frequently  used 
to  compute  the  model  function  and  derivatives  for  each  row,  as  in  the  case 
of  the  example  shown  in  Appendix  A. 

When  the  value  of  the  user-supplied  derivative  disagrees  with  the  cor¬ 
responding  finite  difference  value,  the  checking  procedure  attempts  to  de¬ 
termine  whether  the  disagreement  is  due  to  am  error  in  the  user’s  code  or 
is  due  to  the  inaccuracy  of  the  finite  difference  approximation.  The  check¬ 
ing  procedure  automatically  generates  an  error  report  when  one  or  more 
of  the  derivatives  are  found  to  be  either  incorrect  or  questionable.  This 
report  identifies  which  derivatives  appear  correct,  which  appear  incorrect, 
and  which  appear  questionable  and  for  what  reason. 
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2.2.4  Scaling 

Poorly  scaled  problems,  i.e.,  problems  in  which  the  unknowns,  0  and  6,  vary 
over  several  orders  of  magnitude,  can  cause  difficulty  for  least  squares  proce¬ 
dures.  ODRPACK’s  scaling  algorithm  attempts  to  automatically  overcome 
these  difficulties,  although  it  is  preferable  for  the  user  to  choose  the  units 
of  the  variable  space  so  that  each  of  the  parameters  0  will  have  roughly 
the  same  magnitude  (DenS83].  When  the  variables  have  roughly  the  same 
magnitude,  the  ODRPACK  scaling  algorithm  will  select  scale  values  which 
are  equal,  and  the  resulting  computations  will  be  the  same  (except  for  the 
effect  of  finite  precision  arithmetic)  as  an  unsealed  analysis,  i.e.,  an  analysis 
in  which  all  of  the  scale  values  are  set  to  one.  If  the  variables  do  not  have 
roughly  the  same  magnitude,  the  ODRPACK  scaling  algorithm  will  select 
varying  scale  values.  This  will  not  change  the  optimal  solution  but  may 
affect  the  number  of  iterations  required  to  find  the  solution,  and,  in  some 
cases,  whether  the  algorithm  is  or  is  not  successful. 

Users  may  also  select  their  own  scale  values  as  discussed  in  the  ODR¬ 
PACK  reference  guide. 

2.2.5  Default  Values  and  Structured  Arguments 

ODRPACK  uses  default  values  and  structured  arguments  to  simplify  the 
user  interface.  The  availability  of  default  values  in  ODRPACK  means  that 
the  user  does  not  have  to  be  concerned  with  determining  values  for  many  of 
the  ODRPACK  arguments  unless  the  problem  being  solved  requires  the  use 
of  non-default  values.  Structured  arguments,  described  below,  can  reduce 
the  amount  of  storage  space  required  for  some  arguments  and  the  work 
required  by  the  user  to  initialize  those  arguments. 

Default  Values.  Default  values  have  been  specified  for  ODRPACK  sub¬ 
program  arguments  wherever  feasible.  These  default  values  are  invoked  by 
setting  the  corresponding  argument  to  any  negative  value.  Arrays  with 
default  values  are  invoked  by  setting  the  first  element  of  the  array  to  a 
negative  value,  in  which  case  only  the  first  value  of  the  array  will  ever  be 
used.  This  allows  a  scalar  to  be  used  to  invoke  the  default  values  of  arrays, 
thus  saving  space  and  eliminating  the  need  to  declare  such  arrays. 


6 


Users  are  encouraged  to  invoke  the  default  values  of  arguments  wherever 
possible.  These  values  have  been  found  to  be  reasonable  for  a  wide  class  of 
problems.  Their  use  will  greatly  simplify  the  initial  use  of  ODRPACK  for 
a  given  problem.  Fine  tuning  of  these  arguments  can  then  be  done  later  if 
necessary. 

Structured  Arguments.  Certain  ODRPACK  arguments  specify  attrib¬ 
utes  of  the  individual  components  of  X  and  A,  where  X  is  the  n  x  m  array 
with  t-th  row  x,  and  A  is  the  corresponding  n  x  m  array  with  t-th  row  6,. 
These  attribute  arrays  are  frequently  either  constant  for  all  components  or 
are  constant  within  each  column  and  vary  only  among  the  columns.  The 
structured  argument  facility  allows  a  scalar  to  specify  an  attribute  of  sun 
entire  column  or  of  the  whole  array. 

For  example,  one  such  attribute  array  which  might  display  this  structure 
is  the  matrix  of  6  weights,  RHO ,  where  RHO  is  an  n  x  m  array  with  i-th 
row  the  diagonal  elements  of  A  (c f.  (2)).  Suppose  each  row  of  X  indicates 
an  hourly  temperature  reading  and  each  column  a  different  day  on  which 
the  temperature  readings  were  taken.  Then  the  user  would  probably  want 
to  weight  each  component  of  A  equally,  and  thus  RHO  would  be  constant 
throughout.  If  one  column  of  the  independent  variable  contained  hourly  . 
temperature  readings  and  the  other  hourly  humidity  readings,  then  the 
user  might  want  to  weight  each  component  in  the  first  column  of  A  the 
same,  and  to  weight  each  component  in  the  second  column  the  same,  but 
not  necessarily  want  to  weight  the  two  columns  equally.  In  this  case,  the 
components  of  RHO  would  be  constant  within  each  column  and  would  vary 
only  among  the  columns.  Of  course,  in  other  cases,  the  user  might  want 
to  weight  each  component  of  A  differently,  and  each  component  of  RHO 
would  be  different. 

ODRPACK  structured  arguments  exploit  this  structure.  If  each  of  the 
n  x  m  elements  of  an  attribute  array  is  the  same,  then  a  single  value  can 
be  used  to  specify  all  n  x  m  elements.  If  the  values  of  such  an  array  only 
vary  among  the  columns,  then  each  column  of  the  array  can  be  specified 
by  a  single  value.  Thus,  it  is  only  necessary  to  supply  all  n  x  m  elements 
of  the  structured  argument  array  when  the  elements  of  one  or  more  of  the 
columns  must  be  individually  specified. 
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3.2.6  Calling  Sequences 

As  we  noted  in  §  2.1,  ODRPACK  has  two  levels  of  user 'control  of  the  com¬ 
putations.  These  levels  are  provided  by  user-callable  subprograms  SODR 
and  SODRC  in  single  precision,  and  DODR  and  DODRC  in  double  pre¬ 
cision.  SODR  and  DODR  have  shorter  call  statements  than  SODRC  and 
DODRC  because  they  preset  many  variables  to  their  default  values  in  order 
to  reduce  user  input.  (See  §  2.2.5.)  We  believe  that  SODR  and  DODR  are 
adequate  for  most  problems.  SODRC  and  DODRC,  with  their  expanded 
call  statements,  provide  the  user  with  maximum  flexibility  in  finding  the 
weighted  orthogonal  distance  solution. 

2.2.7  Automatic  Output 

ODRPACK  automatically  generates  computation  reports,  thus  saving  the 
user  from  the  tedious  chore  of  formatting  the  results  for  output.  The 
ODRPACK  computation  reports  are  divided  into  three  parts:  an  initial 
summary,  iteration  summaries  and  a  final  summary.  Each  part  can  be 
generated  in  either  a  long  or  short  form,  and  the  frequency  of  the  itera¬ 
tion  summaries  can  be  specified.  By  default,  the  computation  reports  in¬ 
clude  a  “long”  initial  summary,  no  iteration  summaries  and  a  “short”  final 
summary,  as  shown  in  Appendix  B.  ODRPACK  user-callable  subprograms 
SODR  and  DODR  automatically  generate  the  default  report;  subprograms 
SODRC  and  DODRC  allow  the  user  to  control  the  computation  reports. 

2.3  Example  Program 

Appendix  A  shows  an  example  of  a  user-supplied  driver  for  invoking  ODR¬ 
PACK.  This  sample  program  uses  the  ODRPACK  double  precision  sub¬ 
program  DODR  to  solve  exercise  I  on  pages  521  -  522  of  (DraS81j.  (An 
example  using  ODRPACK  single  precision  subprogram  SODR  could  easily 
be  generated  from  the  one  shown  in  Appendix  A  by  substituting  SODR  for 
DODR  and  changing  the  DOUBLE  PRECISION  declaration  statement  to 
REAL.)  The  model  for  this  example  is 

* =  “p  | _A(*U  +  <,i)Ktp  [-*  -  as)] )  - 
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W. 


The  report  generated  by  ODRPACK  for  this  example  is  shown  in  Appendix 
B. 

3  Installing  ODRPACK 

3.1  Supplied  Software 

3.1.1  ODRPACK  Code 

ODRPACK  software  is  available  in  both  single  and  double  precision  ver¬ 
sions.  These  are  distinguished  by  the  first  letter  of  the  subprogram  name: 

D  for  double  precision,  S  for  single  precision.  The  code  for  each  version  is 
separated  into  three  sections  to  facilitate  installation. 

The  first  section  includes  all  subprograms  written  especially  for  ODR¬ 
PACK.  The  two  user-callable  ODRPACK  subprograms  of  each  version  are 
listed  first,  followed  by  the  ODRPACK  exercise  subprogram.  The  remain¬ 
ing  subprograms  are  then  listed  in  alphabetical  order.  The  code  in  this 
section  should  not  require  any  modification,  unless  the  installer  wishes  to 
customize  the  user-callable  subprograms  as  discussed  in  §  3.2. 

The  second  section  of  code  includes  the  subprograms  used  by  ODR¬ 
PACK  from  the  public  domain  packages  LINPACK  (DonMBS79]  and  BLAS  * 
[LawHKK79j,  also  listed  in  alphabetical  order.  The  installer  can  use  local 
versions  of  these  packages  if  available.  This  would  be  particularly  beneficial 
if  the  installer’s  machine  has  specially  optimized  versions  of  LINPACK  or 
BLAS. 

The  third  section  includes  the  only  machine  dependent  subprogram  in 
ODRPACK.  This  subprogram  supplies  the  machine  dependent  constants 
used  by  ODRPACK.  It  already  includes  comment  statements  listing  the 
necessary  constants  for  a  number  of  common  machines.  If  the  constants 
for  the  target  machine  are  included  in  this  subprogram,  then  the  installer 
need  only  “uncomment”  the  appropriate  DATA  statements.  Note  that  this 
subprogram  will  return  an  undefined  value  until  it  is  updated;  the  installer 
must  update  it  before  compiling  and  running  ODRPACK. 


3.1.2  Test  Drivers  and  Data 


The  ODRPACK  supplied  software  includes  driven  and  data  sets  for  running 
ODRPACK  in  both  single  and  double  precision.  There  are  three  driven 
for  each  venion  of  the  code.  Two  of  these  driven  are  simple  programs 
that  users  can  modify  to  form  their  own  ODRPACK  driven.  The  third  is 
a  demonstration  program  which  exercises  ODRPACK 's  main  features  and 
can  be  used  to  verify  that  the  installation  was  completed  successfully. 

The  demonstration  program  calls  ODRPACK  several  times,  with  each 
call  testing  one  or  more  features  of  the  package.  The  results  of  each  call  are 
automatically  compared  to  the  results  obtained  by  the  authon  using  the 
double  precision  venion  of  ODRPACK  run  on  a  CYBER  855  under  NOS 

2.4.2  (120  bits  per  double  precision  value).  The  success,  or  failure,  of  each 
test  is  noted  individually  in  the  output  generated  by  this  program,  and  is 
summarized  for  all  of  the  data  sets.  By  running  this  demonstration,  the 
installer  can  easily  be  assured  that  the  package  is  performing  as  it  should. 

3.2  Customizing  ODRPACK 

ODRPACK  is  designed  to  easily  accommodate  a  wide  variety  of  problems 
without  any  modification  to  the  supplied  code.  However,  installen  might  . 
want  to  customize  the  user  interface,  i.e.,  the  user-callable  subprograms,  or 
to  customize  the  ODRPACK  generated  reports  for  specific  tasks. 

Customizing  ODRPACK  is  simplified  by  the  extensive  comments  con¬ 
tained  within  each  subprogram,  and  the  consistent  use  of  variable  names 
between  subprograms.  The  SLATEC  Source  File  Format  |FonJS82]  has 
been  followed,  and  each  subprogram  provides  a  standardized  prologue  de¬ 
scribing  the  purpose  of  the  subprogram  and  what  other  subprograms  are 
called,  an  alphabetical  list  of  all  variables  referenced  by  the  subprogram 
and  how  they  are  used,  as  well  as  comments  explaining  the  major  sections 
of  the  code. 

In  addition,  ODRPACK’s  automatically  generated  reports  are  produced 
by  separate  subprograms.  Because  the  report  generators  are  isolated  in  this 
manner,  it  is  relatively  easy  to  produce  customized  reports  if  required. 


3.3  Portability 

ODEPACK  code  conforms  to  the  ANSI  X3  9-1978  FORTRAN  itandird 
and  has  ban  succ— fully  installed  on  a  Cyber  665,  Cyber  306,  IBM  PC/AT, 
Concurrent  3330  and  Vax  11/710.  We  believe  it  will  be  passible  to  install 
ODEPACK  on  any  system  with  an  ANSI  Fortran  77  compiler  and  adequate 


A  Sample  program 

PROGRAM  8 AMPLE 

C 

C  SET  PARAMETERS  FOR  MAXIMUM  PROBLEM  SIZE  HANDLED  BY  THIS  DRIVER 
C  WHERE  MAXI  IS  THE  MAXIMUM  HUMBER  OF  OBSERVATIONS  ALLOWED. 

C  MAXM  IS  THE  MAXIMUM  NUMBER  OF  COLUMNS  IN  THE 

C  INDEPENDENT  VARIABLE  ALLOWED.  AND 

C  MAXNP  IS  THE  NAXINUN  NUMBER  OF  FUNCTION  PARAMETERS 

C  ALLOWED. 

C 

INTEGER  MAXN, MAXM. MAXNP 
PARAMETER  (NAXM-16 .NAXM-6.MAXMP-6) 

C 

C  SET  PARAMETERS  TO  SPECIFY  DIMENSIONS  OF  ARRAYS  USED  BY  DODR 
C 

INTEGER  LDX.LDRHO.LWORK.LIWORX 
PARAMETER 

♦  (LDX-MAIM, 

♦  LDRHOM. 

♦  LWORK  ■  9  ♦  ••MAXM  ♦  10*NAXM*NAXN  ♦  2* MAIN* MAXNP  ♦  8*NAXNP. 

♦  LI WORK  -  2* MAXNP  ♦  MAXM  *  10) 

C 

C  DECLARE  U8ER-SUPPLIED  SUBROUTINES  AND 
C  ALL  OTHER  NECESSARY  VARIABLES  AND  ARRAYS 

C 

EXTERNAL 

♦  FUN. JAC 
INTEGER 

♦  N.N.MP. 

♦  JOB. 

♦  IWORK(LIWORK). 

♦  INFO 

DOUBLE  PRECISION 

♦  X(LDX.NAXM). 

♦  Y(LDX) . 

♦  BETA (MAXNP) . 

♦  RMO(LDRNO.MAXM) . 

♦  WORX(LWORX) 


c 

OPEN (UNIT-5 .FILE- 'DATA1 ' ) 

OPEN (UNIT-C .FILE- 'REPORT ' ) 

C 

C  READ  NUMBER  OF  OBSERVATIONS 

C  NUMBER  OF  COLUMNS  OF  DATA  IN  THE  INDEPENDENT  VARIABLE 

C  NUMBER  OF  PARAMETERS 

C  OBSERVED  VALUES  OF  INDEPENDENT  AND  DEPENDENT  VARIABLES 

C  STARTING  VALUES  OF  FUNCTION  PARAMETERS 

C 

READ  (5.*)  N.M.NP 

READ  (5.*)  ((X(I.J).I-l.N).J-l.M) 

READ  (6.*)  (T(I).I-l.N) 

READ  (6.*)  (BETA(l).I-l.NP) 

C 

C  SPECIFY  DELTA  WEIGHTS 
C 

RHO (1.1)  -  3.0D0 
RH0(1.2)  •  5.0D0 
C 

C  8ET  JOB  TO  DEFAULT  SETTING.  INDICATING 
C  -  SOLUTION  TO  BE  FOUND  BY  ODR 

C  -  DERIVATIVES  TO  BE  COMPUTED  BY  FINITE  DIFFERENCES 

C  -  DELTA'S  TO  BE  INITIALIZED  TO  ZERO 

C 

JOB  -  -1 
C 

C  COMPUTE  ODR  80LUTI0M  USING  FINITE -DIFFERENCE  DERIVATIVES 
C 

CALL  DOOR 

♦  (FUN,  JAC. 

♦  N.M.NP. 

♦  X.LDX, 

♦  T. 

♦  BETA . 

♦  RHO , LDRHO , 

♦  JOB. 

♦  WORE . LWORX . I WORE . LIVORX , 

♦  INFO) 
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END 

SUBROUTINE  FUN(N,NP ,N .BETA . XPLUSD . LDXPD.F, IFU6) 

C 

C  INPUT  ARGUMENTS 

C  (WHICH  MUST  NOT  BE  CHANGED  BY  THIS  ROUTINE) 

C 

INTEGER  N.NP.N.LDXPD 

DOUBLE  PRECISION  BETA(NP) . XPLUSD (LDXPD.M) 

C 

C  OUTPUT  ARGUMENTS 
C 

DOUBLE  PRECISION  F(N) 

INTEGER  I FLAG 

C 

DO  10  I  »  1.  N 

F(I)  •  EXP(-BETA(1)*XPLUSD(I , 1)* 

♦  EXP(-BETA(2)* 

♦  (1.0D0/XPLU8DCI.2)  -  1 . 0D0/620 . ODO) ) ) 
10  CONTINUE 

C 

RETURN 

END 


B  Sample  output 


Report  generated  by  DODR  example  program  from  Appendix  A,  run  using 
a  Cyber  160/855  under  NOS  2.4.2. 


•  ODRPACK  VERSION  1.2  OF  12-16-66  (DOUBLE  PRECISION)  * 


INITIAL  SUMMARY  FOR  FIT  BY  METHOD  OF  OCR 


PROBLEM  SIZE: 


NUMBER  OF  OBSERVATIONS 

NUMBER  OF  COLUMNS  OF  DATA  IB  INDEPENDENT  VARIABLE 
NUMBER  OF  FUNCTION  PARAMETERS 
NUMBER  OF  UNFIXED  FUNCTION  PARAMETERS 


8 

2 

2 

2 


£ 


INDEPENDENT  VARIABLE  AND  DELTA  WEIGHT  SUMMARY: 


X  - 
FIXED  - 
INITIAL  DELTA  • 


8 

i 

a* 

COLUMN  2 

3 

OSS  1 

OBS  N 

OBS  1 

OBS  N 

£ 

. 106000+03 

.660000*02 

.600000*03 

.640000*03 

V 

NO 

NO 

NO 

NO 

a' 

.  000000*00 

.000000*00 

15 

.000000*00 

.000000*00 

9 

vv 

i 

DELTA  SCALE  *  .787400-03  787400-03  .168360-02  .166360-02 

DELTA  WEIGHTS  -  .300000*01  300000*01  600000*01  .600000*01 


OKPnonrr  vaxiailx  amd  obsebvatiomal  mot  weight  sunmaby: 


088  1  088  8 
T  -  .813000*00  376000*00 

088.  18108  VT8.  -  .100000*01  100000*01 


FUNCTION  PAEAMETEB  80MU8T : 


I8DCX  -  1  2 

18ITIAL  8BTA  -  116600000-01  600000000*04 

riXID  -  80  80 

BETA  SCALK  -  .866800870*03  30000000003 


C08T80L  VALUK8  A8D  f TOP PI  IK  CtITCAIA 


JOB  TAUT AC  88T0L  PAATOL  MAX IT 

-0001  100*01  .600-14  860-10  60 

• 

A.  FIT  18  BOT  A  KCSTAKT . 

B  DELTAS  ABE  I8ITIALIZED  TO  ZEBO 
C.  DEEIVAT1VE8  ABE  COMPUTED  BY  FIBITE  DIFFEBE8CE8 . 

0  FIT  IS  BY  HETNOO  OF  08T8000MAL  DISTANCE  B1CBE8SI08 
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INITIAL  SUNS  OP  SQUARES: 


.676620110*00 
.  000000000*00 
.676620110*00 


SUN  OP  SQUAISD  NIIGHTKD  OBSERVATIONAL  ERRORS 

SUN  OP  SQUARED  WEIGHTED  DELTAS 

SUN  OP  SQUARED  VEI68TED  EPSILONS 


PINAL  SUNNART  POR  PIT  BT  NETHOO  OP  ODR 


8 TOPPING  CONDITION  (INTO  -  1): 


THE  RELATIVE  CHANGE  IN  THE  SUN  OP  THE  SQUARED 
WEIGHTED  OBSERVATIONAL  ERRORS  IS  LESS  THAN  SSTOL 

CONDITION 

NUNBER  OP  NUNBER  OP  RUBBER  RANK 

ITERATIONS  PN  EVAL8  (INVERSE)  DEPICIENCT 
6  M  . 1S99D-06  0 


PINAL  SUM  OP  SQUARES: 


SUN  OP  SQUARED  WEIGHTED  OBSERVATIONAL  ERRORS  .763623230-03 
SUN  OP  SQUARED  WEIGHTED  DELTAS  . 236420090-07 
SUN  OP  SQUARED  WEIGHTED  EPSIL0N8  .763796690-03 
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ESTIMATED  BETA(I).  I  *  1 . NP: 


INDEX  VALUE . > 

1  TO  2  .306707270-02  .276273270*06 


ESTIMATED  EPSILON(I)  AND  DELTA(I ,*) .  1*1 . N: 


I  EPSILON(I)  DELTA (I .1)  DELTA(I .2) 

1  .167624450-02  .140661720-06  .424167060-06 

2  .204347180-02  .128362220-05  .202628100-06 

3  -.206000660-01  -.716622010-06  -.233688240-04 

4  .243068320-02  .160470920-05  .241144810-06 

6  .727774820-02  .233032610-06  .620703130-06 

6  .407032640-02  .241628460-06  .404836520-06 

7  .130430710-01  .433373440-06  .147267270-04 


8  -.864006400-02  -.613046600-06  -.848610630-06 
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